#------------------------------------------------------------------------------
rm(list=ls())
library(LalRUtils)
libreq(tidyverse, data.table, zoo, tictoc, fst, 
       fixest, PanelMatch, patchwork,
       rio, magrittr, janitor, did, panelView, 
       ggiplot, tictoc, binsreg, interflex)
set.seed(42)
theme_set(lal_plot_theme())

#------------------------------------------------------------------------------





#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])
#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------

# %%
if (exists("vcf_data")){
  vcf = vcf_data[year>= 1995]
}
vcf[, t := year - 1995]
vcf[, time := year - first_pesa_exposure]

(ex_ante_med = quantile(vcf$cover_1990, 0.5))
above_med = vcf[cover_1990 > ex_ante_med]
above_med[, never_treated := max(D) == 0, cellid]


# %% jharkhand - VCF
above_med[, jh := ifelse(state == "Jharkhand", "Jharkhand", "Others")]
bystate_vcf = feols(forest_index ~ D | cellid[t] + styear,
                    data = above_med, cluster = "blk", fsplit = ~ jh)

controls_pre1 = above_med[never_treated == 1 & year < first_pesa_exposure, 
                          mean(forest_index)]
controls_pre2 = above_med[never_treated == 1 & state == "Jharkhand" & 
                            year < first_pesa_exposure, 
                          mean(forest_index)]
controls_pre3 = above_med[never_treated == 1 & jh == "Others" & 
                            year < first_pesa_exposure, 
                          mean(forest_index)]

treat_pre1    = above_med[never_treated == 0 & 
                            year < first_pesa_exposure, 
                          mean(forest_index)]
treat_pre2    = above_med[never_treated == 0 & state == "Jharkhand" & 
                            year < first_pesa_exposure, 
                          mean(forest_index)]
treat_pre3    = above_med[never_treated == 0 &  jh == "Others" & 
                            year < first_pesa_exposure, 
                          mean(forest_index)]


# %%
gfc[, jh := ifelse(state == "Jharkhand", "Jharkhand", "Others")]
bystate_gfc = feols(def_ha ~ D | village + village[t] + styear, 
                    cluster = "block",
                    gfc[gfc3 == TRUE & pref == 1], fsplit = ~ jh)

gfc[, ever_treated := max(D), .(village)]
ctrl_mean1  = round( gfc[ pref == 1 & gfc3 == T & ever_treated == 0 & 
                         pesa_exposure == 0, mean( def_ha ) ], 2 )
ctrl_mean2  = round( gfc[ pref == 1 & gfc3 == T & ever_treated == 0 & 
                            state == "Jharkhand" & 
                           pesa_exposure == 0, mean( def_ha ) ], 2 )
ctrl_mean3  = round( gfc[ pref == 1 & gfc3 == T & ever_treated == 0 & 
                            jh == "Others" & 
                           pesa_exposure == 0, mean( def_ha ) ], 2 )

treat_mean1 = round(gfc[pref == 1 & gfc3 == T & ever_treated == 1 & 
                          pesa_exposure == 0, 
                        mean(def_ha)], 2)
treat_mean2 = round(gfc[pref == 1 & gfc3 == T & ever_treated == 1 & 
                          state == "Jharkhand" & pesa_exposure == 0, 
                        mean(def_ha)], 2)
treat_mean3 = round(gfc[pref == 1 & gfc3 == T & ever_treated == 1 & 
                          jh == "Others" & pesa_exposure == 0, 
                        mean(def_ha)], 2)
# %%
treatmap =c(
  # GFC
  "def_ha"       = "Annual Deforestation in Hectares",
  "village"      = "Village",
  "village[t]"   = "Village + Village TT",
  # VCF
  "forest_index" = "Forest cover index",
  "green_index"  = "Non-forest green index",
  "built_index"  = "Non-forest index",
  "cellid"       = "Pixel",
  "cellid[t]"    = "pixel + pixel TT",
  # both
  "yr"           = "Year",
  "year"         = "Year",
  "styear"       = "State $\\times$ Year",
  "D"            = "PESA $\\times$ Scheduled",
  "block"        = "Block",
  "blk"          = "Block"
)

# %%
etable(bystate_vcf, bystate_gfc)


# %%
fitstat_register("n_new", function(x) summary( x )$nobs , 
                 "\\# Observations")
fitstat_register("clust", function(x) "Block" , 
                 "Standard-Errors")
etable(bystate_vcf, bystate_gfc,
       style.tex = style.tex(main = "base", 
                             depvar.title = "", 
                             model.title = "", 
                             var.title = "\\midrule", 
                             slopes.title = "\\midrule \\emph{Time Trends}", 
                             yesNo = c("$\\checkmark$", ""), 
                             signif.code = NA, 
                             line.bottom = "\\midrule \\midrule"),
       headers = c( "Full sample", "Jharkhand", "Others", 
                    "Full sample", "Jharkhand", "Others" ), 
       signif.code = NA,
       fixef_sizes = T, 
       fixef_sizes.simplify = F,
       fitstat = ~ n_new + r2, 
       tex = TRUE,
       dict = treatmap,
       label = "tab:regs_all_bystate",
       title = glue::glue("Regression estimates decomposed by 
                          state (analysis sample at ex-ante median cutoff)"),
       extralines=list(
         "\\midrule \\emph{Summary Statistics}" = c("", "","","","","" ),
         "Mean Pre-Y (Non-Sch)"= c( controls_pre1, controls_pre2, 
                                    controls_pre3, ctrl_mean1, 
                                    ctrl_mean2, ctrl_mean3 ),
         "Mean Pre-Y (Sch )"   = c( treat_pre1, treat_pre2, treat_pre3, 
                                treat_mean1, treat_mean2, treat_mean3 ),
         "Dataset"     = c( rep( "VCF", 3 ), rep( "GFC", 3 ) ),
         "Timespan"    = c( rep( "1995-2017", 3 ), rep( "2001-2017", 3 ) )
       ),
       file = "appendix_tableA5.tex", 
       replace = TRUE
)

#------------------------------------------------------------------------------
